Phase space structure and dynamics for the Hamilt 

isokinetic thermostat 



Peter Collins/ Gregory S. Ezra,^'!^] and Stephen Wiggins^'j^ 

^School of Mathematics 
University of Bristol 
Bristol BS8 ITW 
United Kingdom 
^Department of Chemistry and Chemical Biology 
Baker Laboratory 
Cornell University 
Ithaca, NY 14853 
USA 

(Dated: April 27, 2010) 



1 



Abstract 

We investigate the phase space structure and dynamics of a Hamiltonian isokinetic thermostat, 
for which ergodic thermostat trajectories at fixed (zero) energy generate a canonical distribution 
in configuration space. Model potentials studied consist of a single bistable mode plus transverse 
harmonic modes. Interpreting the bistable mode as a reaction (isomerization) coordinate, we es- 
tablish connections with the theory of unimolecular reaction rates, in particular the formulation of 
isomerization rates in terms of gap times. The distribution of gap times (or associated lifetimes) 
for a microcanonical ensemble initiated on the dividing surface is of great dynamical significance; 
an exponential lifetime distribution is usually taken to be an indicator of 'statistical' behavior. 
Moreover, comparison of the magnitude of the phase space volume swept out by reactive trajecto- 
ries as they pass through the reactant region with the total phase space volume (classical density 
of states) for the reactant region provides a necessary condition for ergodic dynamics. We compute 
gap times, associated lifetime distributions, mean gap times, reactive fluxes, reactive volumes and 
total reactant phase space volumes for model systems with 3 degrees of freedom, at three different 
temperatures. At all three temperatures, the necessary condition for ergodicity is approximately 
satisfied. At high temperatures a non-exponential lifetime distribution is found, while at low tem- 
peratures the lifetime is more nearly exponential. The degree of exponentiality of the lifetime 
distribution is quantified by computing the information entropy deficit with respect to pure ex- 
ponential decay. The efficacy of the Hamiltonian isokinetic thermostat is examined by computing 
coordinate distributions averaged over single long trajectories initiated on the dividing surface. 
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I. INTRODUCTION 



Deterministic thermostats are widely used to simulate equilibrium physical systems de- 
scribed by ensembles other than microcanonical (constant energy and volume, {E, V)), such 
as constant temperature- volume iT,V) or temperature-pressure (T,p) [IHE], or for simula- 
tions of nonequilibrium phenomena [7HT2]. 

Deterministic thermostats are obtained by augmenting the phase space variables of the 
physical system of interest with a set of additional variables whose role is to alter the 
standard Hamiltonian system dynamics in such a way that a suitable invariant measure 
in the system phase space is preserved [H [T31 [E]. For example, in the familiar Nose- 
Hoover (NH) thermostat flSl [H] the exact dynamics preserves both an extended energy 
and a suitable invariant measure, ensuring that, provided the extended system dynamics is 
effectively ergodic on the timescale of the simulation, the physical system will sample its 
phase space according to the canonical (constant T) measure. 

Extended system thermostat dynamics can be either Hamiltonian [T3HTH] or non- 
Hamiltonian [TJ [THHSS]- For example, NH dynamics is non-Hamiltonian [T3]. (The NH 
system is however conformally symplectic [271 [28] , meaning that the dynamics can be ren- 
dered Hamitonian by a coordinate-dependent scaling of the vector field |27] . This is not true 
for other non-Hamiltonian thermostats such as Nose- Hoover chains [29] or Bulgac-Kusnezov 
global demons [50] .) 

There has been considerable interest in the formulation of Hamiltonian deterministic ther- 
mostats such as the Nose-Poincare system [18] . In this approach, the extended Hamiltonian 
for the physical system plus thermostat variables incorporates a coordinate-dependent time 
scaling of Poincare-Sundman type [311 ES] . Restricting the dynamics to a fixed value (zero) 
of the extended Hamiltonian results in the system variables sampling their phase space ac- 
cording to (for example) the canonical density [IB] (again, subject to the assumption of 
ergodicity). An important motivation for the introduction of Hamiltonian thermostats is 
the possibihty of using symplectic integration algorithms to integrate trajectories [Hl33]IM]. 

As already indicated, a fundamental question concerning deterministic thermostats has 
to do with the effective ergodicity of the dynamics on the timescale of the simulation. 
If the dynamics is not effectively ergodic, then trajectory simulations will not generate 
the correct invariant measure [3S1 ES]- It has long been recognized, for example, that the 
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dynamical system consisting of a single harmonic oscillator degree of freedom coupled to 
the NH thermostat variable is not ergodic [E] (see also refs IH7H5^ : in fact, for typical 
parameter values the system phase space exhibits extensive persistence of invariant tori 
(quasiperiodic motion) jlOl SI] . A large amount of effort has been expended in attempts to 
design thermostats exhibiting dynamics more ergodic than the basic NH system [U [5l [29], Il2] . 
(For a careful discussion of the question of ergodicity for 'typical' interatomic potentials, see 
|43j : for related fundamental questions in molecular dynamics, see also ref. 

The question of ergodicity in thermostats is conceptually closely related to the problem 
of statistical versus nonstatistical behavior in the (classical) theory of unimolecular reaction 
rates jlHHlH] • Broadly speaking, in this case one would like to know whether a molecule will 
behave according to a statistical model such as RRKM theory, or whether it will exhibit 
significant deviations from such a theory, ascribable to nonstatistical dynamics |l9H5l] . Such 
'nonstatisticality', which can arise from a number of dynamical effects, can be considered to 
be analogous to the failure of ergodicity in deterministic thermostats. 

In recent years there have been significant theoretical and computational advances in 
the application of dynamical systems theory [551 - I57] to study reaction dynamics and phase 
space structure in multimode models of molecular systems, and to probe the dynamical ori- 
gins of nonstatistical behavior [58] - [63] . The fundamental chemical concept of the transition 
state, defined as a surface of no return in phase space, has been successfully and rigorously 
generalized from the well-established 2 degree of freedom case [M] to systems with > 3 
degrees of freedom Moreover, dynamical indicators exist (determination of reactive 

phase space volume, behavior of the reactive flux) to diagnose nonstatistical behavior (see, 
for example, Ref. [SH [SI [65] ) . 

Nevertheless, relatively little work has been done applying the powerful techniques from 
dynamical system theory [551457] to study the phase space structure and dynamics of deter- 
ministic thermostats. Hoover and coworkers have investigated the fractal nature of various 
phase space structures for equilibrium [66] and nonequilibrium [67I - I69] versions of the NH 
thermostat. Using a Hamiltonian formulation of the isokinetic thermostat [I5] (cf. Section 
[H] of this paper), Morriss and Dettmann [2] have mapped the dynamics of the isokineti- 
cally thermostatted Lorentz gas onto the geodesic motion of a free particle on a particular 
Riemannian manifold. Leimkuhler and Sweet have used Poincare surfaces of section and 
dynamical frequency analysis in their analysis of optimal coupling constants for NH and 
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related thermostats jl2] , while Legoll and coworkers have applied KAM theory to rigorously 
prove the existence of invariant tori in the NH thermostatted harmonic oscillator [101 HI] . 
D'Alessandro, Tenenbaum and Amadei have computed Lyapunov exponents for a ther- 
mostatted united-atom model of the butane molecule [TOj; both Nose- Hoover and isokinetic 
Gaussian thermostats were considered. 

Thermostatted systems present several difficulties for detailed studies from a phase space 
perspective. First, such systems usually have high dimensionality; for example, a NH chain 
thermostatted system has at least 3 degrees of freedom. For such a system the Poincare 
surface of section analysis, standard for 2 degrees of freedom, does not afford the advantage of 
dimensional reduction and global visualization since the surface of section is four dimensional 
for three degrees of freedom. Second, there is a lack of readily computable diagnostics that 
can be used to establish ergodicity [36], other than comparison of coordinate distributions 
obtained using the given thermostat with those obtained using other methods. 

In the present paper we analyze the Hamiltonian isokinetic thermostat [2] from the per- 
spective of reaction rate theory. Although not as widely used as the NH thermostat and 
its many variants, the non-Hamiltonian version of the isokinetic thermostat has been devel- 
oped and applied to several problems of chemical interest by Minary et al. [711 172] . In this 
thermostat, the particle momenta are subject to a nonholonomic constraint that keeps the 
kinetic energy constant. The resulting dynamics generates a canonical (constant tenpera- 
ture) distribution in configuration space [2]. In this work we investigate a slightly modified 
version of the Hamiltonian isokinetic thermostat given by Litniewski [73] and Morishita [H] . 

The structure of the present paper is as follows: Section [H] reviews the Hamiltonian 
formulation of the isokinetic thermostat. The non-Hamiltonian equations of motion for a 
Hamiltonian system subject to the isokinetic constraint correspond to Hamiltonian dynamics 
at zero energy under an effective Hamiltonian whose potential is obtained from the physical 
potential by exponentiation. The model Hamiltonians analyzed in the present paper are 



introduced in Section |III[ In these systems, the physical potential describes n uncoupled 
oscillators, n — 1 harmonic modes together with a bistable thermalizing [7T1 172] or isomer- 
ization coordinate. For these Hamiltonians the physical potential exhibits a saddle of index 
one, as for the case of a bistable reaction profile coupled to one or more transverse confining 



modes. In Section IV we briefiy review earlier results [73] establishing that the isokinetic 
Hamiltonian dynamical system corresponding to the physical Hamiltonian exhibits a nor- 
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mally hyperbolic invariant manifold (NHIM) associated with the saddle [57], at least for 
a limited range of energies above that of the saddle. The phase space formulation of uni- 
molecular reaction rate theory in terms of the gap time [511 ESHZS] and related concepts is 
discussed in Sect. |Vj The classical spectral theorem [791 - IM] provides a relation between the 
distribution of gap times and the phase space volume occupied by reactive phase points; 
unless the measure of the region swept out by isomer izing trajectories equals that of the 
energy shell, the system cannot be ergodic. This necessary condition for ergodicity estab- 
lishes a connection between the concepts of reaction rate theory and the properties of the 
Hamiltonian isokinetic thermostat. Numerical results on 3 and 4 degree of freedom Hamil- 
tonian isokinetic thermostats are reported in Section VI Section VII| concludes. Analytical 
results for the density of states for exponentiated harmonic potential are given in Appendix 
[K\ while details of the phase space sampling procedures used in our numerical work appear 
in Appendix [B] 

A brief discussion of the theoretical framework developed here was given in ref. [751 
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II. THE HAMILTONIAN ISOKINETIC THERMOSTAT 



The role of the isokinetic thermostat is to dynamically constrain the kinetic energy of 
a simulated system to have a constant value proportional to k-^T, where T is the desired 
temperature and /cb is Boltzmann's constant. Non-Hamiltonian equations of motion for the 
isokinetic thermostat have been obtained by applying Gauss' principle of least constraint 
[21 O [8], which is the appropriate dynamical principle for nonholonomically constrained 
systems [TJ |85] . 

Provided the underlying dynamics is effectively ergodic on the timescale of the simulation, 
a trajectory of the non-Hamiltonian isokinetic thermostat samples configuration space ac- 
cording to the invariant measure associated with the equilibrium canonical ensemble [21 186] . 
Algorithms for the non-Hamiltonan isokinetic thermostat have been developed and applied 
to a variety of systems by Minary, Martyna and Tuckerman [7T| [72] . 

A Hamiltonian formulation of the isokinetic thermostat has been given by Dettmann and 
Morriss [21 US]- The Hamiltonian used in the Dettmann-Morriss approach incorporates a 
coordinate-dependent scaling of time [52] 157] , and the physically relevant dynamics is re- 
stricted to the zero energy hypersurface. A noncanonical tranformation of variables then 
leads to a set of dynamical equations equivalent to the non-Hamiltonian version. An impor- 
tant motivation for the development of Hamiltonian versions of the isokinetic thermostat 
is the possibility of using symplectic integration algorithms to ensure qualitatively correct 
behavior of integrated trajectories over long times [531 El]- 

In this section we briefly review the Hamiltonian formulation of the isokinetic thermostat 
(for more detailed discussion, see [75].) As our focus here is on systems with n >3 degrees of 
freedom (DoF), we discuss a Hamiltonian version of the isokinetic thermostat that generates 
the invariant measure associated with the configurational canonical ensemble, while at the 
same time allowing use of the simplest Verlet-type [88] symplectic integration algorithm 

[ZSlEl]. 

A. Hamiltonian isokinetic thermostat 



The physical Hamiltonian of interest is assumed to have the standard form 
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(1) 



where {q,p) G M" x M", $(g) is the potential energy, and we set all masses equal to unity 
for simplicity. Corresponding physical Hamiltonian equations of motion are 

« = + g^ (2a) 
P^-f. (2b) 

For n ^ 3 degrees of freedom, we define a Hamiltonian Ti = T-L^tt , q) as [73l [H] 

(3) 

where tt G M" and /3 and u are parameters to be determined. The associated Hamiltonian 
equations of motion are 

+ - = . (4a) 
-' = -^ = -.$,6-^ (4b) 

where the time derivative is denoted by a prime (') to indicate that the derivative is actually 
taken with respect to a suitably scaled time variable (see [21 ES])- As the kinetic energy 
term in ^ is coordinate-independent, a basic Verlet-type symplectic integrator can be used 
to integrate Hamilton's equations for "H [U [33l ISl IHH] • (Higher-order symplectic algorithms 
can of course also be used [U IS]-) 

Making a noncanonical change of coordinates from variable (vr, q) to physical variables 
(P,9), (7r,g) (p,g), where 

TT = e-'^^p (5) 



the Hamiltonian function "H becomes 

1 r 7/1 

(6) 
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Setting "H = then automatically enforces an isokinetic constraint in terms of the momen- 
tum variables p [2] 

- = -■ (7) 
2 2/3 ^ ^ 

For trajectories run at = 0, and only for this value, we have the invariant volume 
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element [2] 



dV = d''7rd''q6(n) 



= d>d"g2e- ("-2)^*5 
We therefore set (/3 = I/ZcbT" as usual) 

1 



1^-2/3$ 



/3 



(8a) 
(8b) 

(8c) 
(8d) 

(9) 



{n - 2)k^T {n-2)' 
to ensure that the invariant measure is proportional to the canonical density (recall n > 2) 



p{q) oc exp[-$(g)/A;BT]. 



With 



the kinetic energy 



n 



{n-2) 



~2 



n 



knT 



(10) 



(11) 



(12) 



as required for an n degree of freedom system at temperature T. Nevertheless, the important 
aspect of the dynamics for computing configurational averages is the invariant measure, not 
the magnitude of the constrained KE. The quantity u can then be treated as a free parameter. 



which can be used to move the potential saddle closer to energy £ = (cf. Section IV). In 
fact, we shall take z/ = 1 in our calculations. 

For n > 3 degrees of freedom, and provided the dynamics is ergodic on the H = energy 
surface, Hamiltonian dynamics Q will therefore yield the canonical measure in configuration 
space. Since the aim is to obtain a canonical distribution, the Hamiltonian dynamics derived 
using ([3]) is in this respect equivalent in principle to the original isokinetic thermostat [2]. 

For completeness, we note that, in order to treat the n = 2 DoF case, it is necessary to 
use a different Hamiltonian /C (cf. refs [SlIZS]), with 



/C 



2 2/3 

To obtain the correct canonical measure case we must then take (n ^ 2) 

1 ^ 



(n - 1)A;bT (n - 1) 



(13) 



(14) 



and 



n 



(n-1) 



(15) 
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III. MODEL HAMILTONIANS 



In this section we introduce the model Hamiltonians studied in the remainder of the 
paper. We define systems with three and four DoF suitable for studying the Hamiltonian 
isokinetic thermostat on the zero energy surface H = 0. 



The phase space structure of these systems is discussed in Section IV, while numerical 
computations of gap times, reactive volumes and coordinate distributions obtained using 
the Hamiltonian isokinetic thermostat are described in Sec. IVII 

The systems considered here consist of n uncoupled oscillators: n — 1 harmonic modes 
plus a single bistable mode. Although the n modes are uncoupled in the physical Hamil- 
tonian H{p,q), exponentiation of the potential ^{q) introduces intermode coupling in the 
isokinetic thermostat Hamiltonian 7i. The bistable mode can be interpreted either as a ther- 
malizing degree of freedom (following Minary, Martyna and Tuckerman (MMT) [TH [72]). 
or as a reactive degree of freedom associated with an isomerization process. Adopting the 
latter perspective, we can apply concepts and methods recently developed for understand- 
ing reaction dynamics (in particular, transition state theory) in phase space [511 EBHSS] to 
investigate the problem of thermostat dynamics. 

A. Double-well potential 

The double well potential for the thermalizing mode is taken to be a temperature- 
independent quartic having the following form: 

xiy) = liy'- (16) 

(Note that the Minary, Martyna & Tuckerman version of the double-well potential is effec- 



tively temperature dependent [7T1[72].) The potential (16) has stationary points at ?/ = (a 



maximum) and y = ±a/q;/2 (minima), with corresponding values and — a^/S. Expanding 
xiy) about the minima y = ±A/a/2 we find that the effective frequency for harmonic motion 
in the vicinity of the minima is 

CO = V2^. (17) 
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B. Model Hamiltonians 



1. Three degrees of freedom 

The 3 DoF potential corresponds to a separable bistable (y) plus harmonic bath modes 
(xi,X2): 

$(xi,X2,i/) = + + 2 ) • (1^) 

Setting m = 1 and choosing potential parameters ui = 1, U2 = \/2 and a = 2 we have 

2 4 

$(Xi,X2,2/) = ^+x2+|--y2. (19) 



The isokinetic thermostat Hamiltonian H is then 

r r -)-[ 

(20) 



The origin is an equilibrium point of saddle-center-center type with energy 
H{0, 0, 0, 0, 0, 0) = — and z/ and /3 are parameters that we can vary. 

As mentioned above, the value of the parameter u determines the (constant) value of the 
physical kinetic energy p'^/2. It also determines the energy of the saddle point with respect 
to the zero of energy, "H = 0. To obtain the correct canonical invariant density, it is only 
necessary that be constant; the parameter u is therefore effectively a free parameter in 
addition to the temperature T. 

For the numerical computations to be discussed below we take a = 2, u = 1. Plots of the 
X2 = slice through the physical potential $ for the 3 DoF system and the corresponding 
exponentiated potential {(3 = 1) are shown in Figure [T} 

2. Four degrees of freedom 

The 4 DoF potential corresponds to a separable bistable (y) plus harmonic bath modes 

muj^x'^ mcOnXn mujlx'i, I / a 9\ / 
$(xi, X2, xs, y) = + + + 2 ) • (21) 

Setting m = 1, cji = 1, Ci;2 = \/2, 1^3 = "\/3 and a = 2 we have 

<^{x,,X2,Xs,y) = ^ + xl + ^ + ^~y^. (22) 
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From eq. (|9]), the parameter (3 = (3/2 for n = 4 DoF, so that the isokinetic thermostat 
Hamiltonian l-L is: 

r2 ^2 ^2 ^2 



2 i i 

TT Vr TT TT 



V 

— — exp 



/3 



-/3 



(23) 



The origin is an equihbrium point of saddle-center-center-center type with energy 
7^(0,0,0,0,0,0,0,0) = — |, and v and (3 are parameters that we can vary. We set the 
parameter v = 1. 



In Section VI we present numerical results for both the 4 DoF and 3 DoF systems corre- 
sponding to three values of the temperature: /3 = 1, 3, and 5. 
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IV. MICROCANONICAL PHASE SPACE STRUCTURE: HAMILTONIAN AND 
NON-HAMILTONIAN ISOKINETIC THERMOSTAT 



In this section, we briefly discuss the phase space structure in the vicinity of the saddle- 
center equilibrium points for Hamiltonians ([T]) and (|3|. 
In previous work, we have established the following [75] : 



(a) If the physical Hamiltonian system defined by eq. ([T]) has an equilibrium point of 
saddle-centre-. . .-centre stability type [57] at the origin then the Hamiltonian system 



defined by an extended Hamiltonian of the form (13) corresponding to the isokinetic 
thermostat has a equilibrium point at the origin of saddle-centre-. . .-centre stability 
type. 

(b) The energy of the equilibrium is such that on the zero energy surface of the extended 
Hamiltonian K, the phase space structures present in the physical Hamiltonian also 
exist in the extended Hamiltonian phase space. 

(c) The phase space structures on the zero energy surface corresponding to the Hamil- 
tonian isokinetic thermostat map to phase space structures in the non-Hamiltonian 
thermostatted system obtained by transforming the Hamiltonian equations of motion 
for evolution (vr, q) under K, to equations of motion for non-canonical variables (p, q) 
under the constraint of zero energy. 



Although these results were previously established for the Hamiltonian /C, eq. (13), they 
also hold for Hamiltonian eq. (|3]) as we describe below. 

More precisely, we have the following: assume that the physical Hamiltonian system ([2]) 
has an equilibrium point at = {q*,p*) = (0,0). The energy of this equilibrium point 

is H{0,0) = $(0). The stability of the equilibrium point is determined by the eigenvalues 
of the derivative of the Hamiltonian vector field (or Hessian of the Hamiltonian function) 
evaluated at the equilibrium point. This is given by the 2n x 2n matrix: 

Hess,^, = , (24) 

\ -$99(0) 0„xn J 

where 0„xn denotes the nxn matrix of zeros and id„xn denotes the nxn identity matrix. 
We require the equilibrium point to be of saddle-centre-. . .-centre stability type. This means 

14 



that the 2n x 2n matrix Hess^ys has eigenvalues ±A, ±iuJi, i = 2, . . . ,n where A and uJi are 
real. 

Eigenvalues 7 of Hess^y^ are obtained by solving the characteristic equation det(HesSgy5 — 
7id2nx2n) = 0. The block structure of the 2n x 2n matrix Hess^ys implies that (cf. Theorem 
3 of [89J) 

det(Hess,y, - 7id2nx2n) = det($gg(0) + 7^id„xn) = (25) 

so that the 2n eigenvalues 7 are given in terms of a, the eigenvalues of the n x n Hessian 
matrix $gg(0) associated with the potential, as follows 



7fc,7fc+„ = ±y/-(Tk, k = l,...,n. (26) 

Therefore, if $(g) has a rank-one saddle at g = 0, so that one eigenvalue is strictly nega- 
tive and the rest are strictly positive, then {q,p) = (0,0) is a saddle-centre-. . .-centre type 
equilibrium point for ^ as described above. 

Next, we consider Hamilton's equations associated with the extended Hamiltonian (|3]), 
eq. Q, corresponding to the isokinetic thermostat. It is easy to verify that {q, vr) = (0, 0) is 
an equilibrium point for Q with energy 7^(0,0) = — ^e~^^**^°^. 

Proceeding as above, we linearize Q about (g, vr) = (0, 0) and compute the eigenvalues 
of the matrix associated with the linearization. These are given by: 

7fc, 7fc+n = iv^e-'^*^'') k = l,...,n. (27) 

In other words, the eigenvalues of the linearization of ^ about (g, tt) = (0, 0) correspond to 
the eigenvalues of the matrix associated with the linearization of ^ about {q,p) = (0,0), 
but with each eigenvalue multiplied by the positive constant y/i'e~^'^^^\ Hence, it follows 
that if the potential of the physical Hamiltonian, $(g), has a rank-one saddle at g = 0, so 
that one eigenvalue is strictly negative and the rest are strictly positive, then (g, vr) = (0, 0) is 
a saddle-centre-. . .-centre type equilibrium point for (|4]). Moreover, if the purely imaginary 
eigenvalues in (26) satisfy a non-resonance condition, then the purely imaginary eigenvalues 



in (27) satisfy the same non-resonance condition |75j . 

The equilibrium point (g, vr) = (0,0) has energy 7/(0,0) = — Tj^e~^^*'^°''. However, we 
are only interested in the dynamics on the 7/ = energy surface. All of the phase space 
structure discussed above exists for a certain range of energies above that of the saddle- 
centre-. . .-centre, and will do so on the "H = surface if — ;j^e~^^**-°'' is close enough to 
zero. 
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Putting together the results above, we have the following [75]: 

Suppose the physical Hamiltonian system (|2| has an equilibrium point of saddle- 
centre-. . .-centre stability type at the origin. Then the Hamiltonian system ^ 
corresponding to the isokinetic thermostat for n > 3 DoF also has an equilibrium 
point of saddle- centre-. . .-centre stability type at the origin. Moreover, the energy 
of this equilibrium point can be chosen so that on the energy surface corresponding 
toTi = there exists a normally hyperbolic invariant manifold (NHIM) /57/ , with 
associated stable and unstable manifolds, a "dividing surface" of no-return and 
minimal flux, and a foliation of the reaction region by n-dimensional invariant 
Lagrangian submanifolds 163]. 

Following the general arguments outlined in [75], it is straightforward to show that the 
phase space structure of Q on "H = exists in the non-Hamiltonian thermostatted system 
in the original physical variables. These general results allow us to conclude that, under the 
noncanonical transformation of variables (vr, g) {p,q), 

• The 2n — 1 dimensional invariant energy surface T-L = corresponding to the Hamil- 
tonian isokinetic thermostat maps to a 2n — 1 dimensional invariant manifold for 
the non-Hamiltonian thermostatted system. 

• The 2n — 3 dimensional NHIM, its 2n — 2 dimensional stable and unstable manifolds, 
the n-dimensional invariant Lagrangian submanifolds, and the 2n — 2 dimensional 
dividing surface map to a2n — 3 dimensional NHIM, its 2n — 2 dimensional stable and 
unstable manifolds, n-dimensional invariant submanifolds, and the 2n — 2 dimensional 
dividing surface in the 2n — l dimensional invariant manifold for the non-Hamiltonian 
thermostatted system. 
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V. PHASE SPACE GEOMETRICAL STRUCTURES, UNIMOLECULAR REAC- 
TION RATES AND THERMOSTAT DYNAMICS: GAP TIMES AND REACTIVE 
VOLUMES 

Our analysis of thermostat dynamics will be carried out in phase space, and our approach 
to probing the dynamics, and especially the question of ergodicity, will be based on the 
general formulation of unimolecular reaction rates based upon that originally given by Thiele 
[77] . In their general form the rate expressions derived by Thiele explicitly invoke the 
existence of a phase space dividing surface separating reactants and products; such surfaces, 
discussed by Wigner [90] (see also refs EU |92]), have only recently become amenable to 
direct computation via the use of normal form approaches [58] - [60| 1631 l83t [HU 1931498] . Our 
discussion of Thiele's reaction rate theory theory follows ref. jSH 



A. Phase space dividing surfaces: definition and properties 

Interpreting the bistable thermalizing mode [Til [72] in the model Hamiltonians of Sec- 



tion III as an isomerization coordinate, we see that the thermostat dynamics associated 
with Hamiltonian ([s]) is equivalent to an isomerization reaction at constant energy "H = 0. 
We therefore consider the rate of unimolecular isomerization, at fixed energy, for a system 
described by a time- independent, n degree-of-freedom (DOF) classical Hamiltonian. 

Points in the 2n-dimensional system phase space Ai = M^" are denoted 2 = (vr, g) G A^. 
The system Hamiltonian is 'H(z), and the {2n — 1) dimensional energy shell at energy E, 
T-ii^z) = E, is denoted T^e C A4. The corresponding microcanonical phase space density is 
6{E — Ti^z)), and the associated density of states for the complete energy shell at energy E 
is 

piE) = [ dz diE-niz)). (28) 
Jm 

In Appendix |A| we provide analytical results for p{E) for n-dimensional systems with 
Hamiltonians of the form ([s]), with $(g) an isotropic harmonic potential and n even. 

The first step in the analysis is to define regions of the energy surface corresponding to 
reactant and product. For the isokinetic thermostat Hamiltonians considered here, there 
is a natural divison of phase space into reactant region (y < 0, say) and product region 
?/ > 0. The dividing surface between reactant and product is determined by symmetry to 
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be the codimension-1 surface y = 0, and we shall be concerned with the evaluation of the 
microcanonical reactive flux across this surface. (For fundamental work on reactive flux 
correlation functions and associated relaxation kinetics, see pUHl03j .) 

For multidimensional systems such as polyatomic molecules {n > 3 DoF), it is in general 
not possible to define or compute a dividing surface with desirable dynamical attributes such 
as the no-recrossing property by working in configuration space alone, and a phase space 
perspective is necessary |58H60l |63l ESI [HIMIH]. 



As discussed in Section IV, Hamilton's equations for the thermostat Hamiltonian Ti have 
an equilibrium of saddle-center-. . .-center stability type. The significance of saddle points of 
this type for Hamilton's equations is that, for a range of energies above that of the saddle, 
the energy surfaces have the bottleneck property in a phase space neighborhood near the 
saddle, i.e., the 2n — 1 dimensional energy surface locally has the geometrical structure of 
the product of a 2n — 2 dimensional sphere and an interval, 5^"^^ x /. In the vicinity of 
the bottleneck, we are able to construct a dividing surface depending on E, DS{E), with 
very desirable properties: For each energy in this range above the saddle, DS(£') locally 
"disconnects" the energy surface into two disjoint pieces with the consequence that the 
only way to pass from one piece of the energy surface to the other is to cross DS(-E'). The 
dividing surface has the geometrical structure of a 2n — 2 dimensional sphere, S*^""^, which 
is divided into two 2n — 2 dimensional hemispheres, denoted DSin(-E) and DSout(-£') that are 
joined at an equator, which is a 2n — 3 dimensional sphere, S*^""^. The hemisphere DSin(-E') 
corresponds to initial conditions of trajectories that enter the reaction region while DSout(-£') 
corresponds to initial conditions of trajectories that exit the reaction region, both by passing 
through the bottleneck in the energy surface. The equator 5*^""^ is an invariant manifold of 
saddle stability type, a so-called normally hyperbolic invariant manifold (NHIM) [57]. The 
NHIM is of great physical significance: it is the actual "saddle" in phase space identified as 
the "activated complex" of reaction rate dynamics [531 1104^ I1U5] . 

In the context of microcanonical rates, it has been shown that DSin(-E) and DSout(-£') 
have the essential no-recrossing property and that the flux across them is minimal [95]. We 
denote the directional flux across these hemispheres by 0in(-E) and 0out(-E), respectively, and 
note that 4>in{E) + (poutiE) = 0. The magnitude of the flux is |0in(-E')| = \(j)ont{E)\ = (t>{E). 
Most significantly, the hemisphere DSin/out(-S) is the correct surface across which to compute 
the "exact" flux into/out of the reaction region. 
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B. Phase space volumes and gap times 



The disjoint regions of phase space corresponding to species A (reactant) and B (product) 
will be denoted Ai^ and Mb, respectively |106j . 

As discussed above, the DS can be rigorously defined to be locally a surface of no return 
(transition state). The microcanonical density of states for reactant species A is 



with a corresponding expression for the density of states Pb{E) for product B for the case 
of compact product energy shell. 

For isokinetic thermostat Hamiltonians H, the = energy surface extends to ±00 
in configuration space. Although the phase space volume N{E) is finite for E 0, the 
corresponding density of states p{E) = dN/dE may diverge as E 0. In Appendix [A| we 
derive analytical expressions for p{E) associated with isotropic harmonic potentials ^{q), 
q G M", with n even. While p{E) diverges as — for n = 2, it is finite in this limit for 
n > 4. Analytical expressions for p{E) are not available for the model potentials studied 
here. 

Provided that the fiow is everywhere transverse to DSjn^ out(-£'); those phase points in 
the reactant region A^a that lie on crossing trajectories |511 1107j (i.e., that will eventually 
react) can be specified uniquely by coordinates {q, tt, ip), where (g, tt) G DSin(-E') is a point on 
DSin(-E), the incoming half of the DS, specified by 2{n — l) coordinates {q, tt), and ip is a. time 
variable. (Dividing surfaces constructed by normal form algorithms are guaranteed to be 
transverse to the vector field, except at the NHIM, where the vector field is tangent [3Sll^ .) 
The point z{q, tt, ip) is reached by propagating the initial condition {q,p) G DSin(-E) forward 
for time ip jSH [77]. As all initial conditions on DSin(£^) (apart from a set of trajectories 
of measure zero lying on stable manifolds) will leave the reactant region in finite time by 
crossing DSout(-£'), for each (g, tt) G DSin(-E') we can define the gap time s = s{q,7r), 
which is the time it takes for the incoming trajectory to traverse the reactant region. That 
is, z{q,7t,ip = s{q,p)) G DSout(-£')- For the phase point z{q,7t,ip), we therefore have 
< < s{q,7z). 

The coordinate transformation z {E, ip, q, tt) is canonical [771 ED l82l 1108] . so that the 
phase space volume element is 




(29) 



d'"2 = dEdx/jda 



(30) 
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with dcr = d" ^gd" ^vr an element of 2n — 2 dimensional area on the DS. 

The magnitude of the flux through dividing surface DS(i?) at energy E is given by 



j 



da 



(31) 



where the element of area da is precisely the restriction to the DS of the appropriate flux 
(2n — 2)-form w*^"^^^ /(n — 1)! corresponding to the Hamiltonian vector field associated with 
H{z) [951 1109milj . The reactant phase space volume occupied by points initiated on the 
dividing surface DSin with energies between E and E + dE is therefore [771 EHHEl] 



dE da dip =dE das (32a) 

ioSinCS) Jo ioSinC-B) 

= dE (j){E) s (32b) 

where the mean gap time s is defined as 

' da s (33) 



1 



HE) 



DSin(E) 



and is a function of energy E. The reactant density of states p^{E) associated with crossing 
trajectories only (those trajectories that enter and exit the reactant region |107j ) is then 

pI{E)=cP{E)-s (34) 



where the superscript C indicates the restriction to crossing trajectories. The result (34) is 
essentially the content of the so-called classical spectral theorem [79HM] . 

If all points in the reactant region of phase space eventually react (that is, all points lie 
on crossing trajectories [CT llU7j ) then 

pI{E)=pp^{E), (35) 

so that the crossing density of states is equal to the full reactant phase space density of 
states. Apart from a set of measure zero, all phase points z G A^a can be classified as either 
trapped (T) or crossing (C) [TD7J. A phase point in the trapped region Ai\ never crosses the 
DS, so that the associated trajectory does not contribute to the reactive flux. Phase points 
in the crossing region M.^ do however eventually cross the dividing surface, and so lie on 
trajectories that contribute to the reactive flux. In general, however, as a consequence of the 
existence of trapped trajectories (either trajectories on invariant trapped n-tori [511 llOTj or 
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trajectories asymptotic to other invariant objects of zero measure), we have the inequahty 
[771 [TU71 [TT^ 

pI{E)<p^{E). (36) 



From the perspective of thermostat dynamics, the equahty (35) is a necessary condition 



for ergodicity. If equahty (35) does not hold, then there is a region of phase space of nonzero 
measure that is trapped on the reactant side of the dividing surface. 

li p^{E) < Pa{E), then it is in principle necessary to introduce corrections to statistical 
estimates of reaction rates \WS\ \107\ I112H115] . Numerical results for p^{E) and p{E) for 
the HCN molecule are discussed in 



C. Gap time and reactant lifetime distributions 

The gap time distribution, V{s; E) is of central interest in unimolecular kinetics [TBI 177] : 
the probability that a phase point on DSin(-E) at energy E has a gap time between s and s+ds 
is equal to V{s; E)ds. An important idealized gap distribution is the random, exponential 
distribution 

V{s;E) = k{E)e-''^''> (37) 

characterized by a single decay constant k (where k depends on energy E), with correspond- 
ing mean gap time s = k~^ . An exponential distribution of gap times is taken to be a 
necessary condition for 'statistical' behavior in unimolecular reactions [7614781 [116j . 

The lifetime (time to cross the dividing surface DSout(-E')) of phase point z{q,p,ip) is 
t = s{q,p) — ip, and the corresponding (normalized) reactant lifetime distribution function 
P(t; E) at energy E is ^MM fTT6HTT9] 

P(t; ^) = _ A Prob(t > t'; E) (38a) 



1 

s 



t'=t 

+ 00 

dsV{s;E) (38b) 



where the fraction of interesting (reactive) phase points having lifetimes between t and t + dt 
is F{t;E)dt. 



Equation (38a) gives the general relation between the lifetime distribution and the frac- 



tion of trajectories having lifetimes greater than a certain value for arbitrary ensembles 



[117H119| . Note that an exponential gap distribution (37) implies that the reactant lifetime 
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distribution F(t; E) is also exponential [76 | [TT f [TT611119] : both gap and lifetime distributions 
for realistic molecular potentials have been of great interest since the earliest days of tra- 
jectory simulations of unimolecular decay, and many examples of non-exponential lifetime 
distributions have been found [651 11141 lll7H123j . 



D. Reaction rates and the inverse gap time 
The quantity 



,HHKM(^) - (39) 

is the statistical (RRKM) microcanonical rate for the forward reaction (A — > B) at energy 
E, the ratio of the magnitude of the flux through DSin(£') to the total reactant density 
of states US HE]. 

Clearly, if px{E) = p2(^), then 



kf^™{E) = - (40) 



the inverse mean gap time. In general, the inverse of the mean gap time is 

1 <P{E) 



s Pa 



^RRKM 



Pa{E) 



(41b) 



> A;™. (41c) 

The inverse gap time can then be interpreted as the statistical unimolecular reaction rate 
corrected for the volume of trapped trajectories in the reactant phase space [THl 11031 11071 

[II21E3]- 

In the next section we discuss our numerical calculations of the following quantities for 
the isokinetic thermostats Hamiltonian at constant energy "H = 0: gap time distributions, 
mean gap time s, reactive flux (f){E = 0), reactive volume s(j) and reactant density of states, 
piE = 0). 
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VI. NUMERICAL COMPUTATIONS FOR ISOKINETIC THERMOSTAT 



A. Computations 

In this Section we discuss the computation of various dynamical quantities for the isoki- 



netic thermostat Hamiltonians defined in Section The potentials for these isokinetic 
thermostat Hamiltonians with 3 and 4 DoF have the form of an exponentiated double well 
plus harmonic modes potential. 

As already noted, for the cases examined here, we can exploit the symmetry of the 
potential functions. Thus, the dividing surface (DS) between 'reactant' and 'product' is 
simply defined to be the symmetry plane ?/ = for the double well potential. For more 
general non-symmetric potentials, the dividing surface in phase space can be computed 
using a normal form expansion [63] . 

1. System parameters 

In addition to the frequencies characterizing the modes transverse to the reaction coor- 
dinate, it is necesssary to choose values for the parameters /3 = l/k-^T, a, and u. 

We adopt the following notation for presentation of our results: the system denoted by 
Yifiav has 3 DoF and parameter values indicated, while the system jpau has 4 DoF. We 



present numerical results for the 3 DoF systems H121, H321, H521 (Hamiltonian eq. (20)) 



and 4 DoF systems J121, J321, J521 (Hamiltonian eq. Q). 

In the units used here the height of the barrier to isomerization in the physical potential 
corresponds to /3 = 2. The temperature values studied here therefore span a range of energy 
scales from well below the barrier height (/3 = 5) to well above (/3 = 1). 



2. Computations 

After choosing a set of parameters, we compute the following quantities (further details 
of computational methodology are given in Appendix [B|) : 

(a) Gap time distribution and mean gap time 
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The distribution of gap times is obtained by starting trajectories on the forward di- 
viding surface and propagating them until they cross the backward dividing surface. 
Initial conditions are obtained by uniform random sampling of the dividing surface 
(see Appendix [b|) . A discretized approximation to the gap time distribution, V{s), is 
obtained by binning the gap times for the trajectory ensemble. The mean gap time 
s is calculated as an unweighted average of computed gap times for the trajectory 
ensemble. 

(b) Lifetime distribution 

The lifetime distribution F{t) is obtained from the (discretized) gap time distribution 



by numerical integration (cf. eq. (38)). The form of the lifetime distibution is of 
interest: in particular, deviations from exponentiality are suggestive of "nonstatistical" 
dynamics. 

The average lifetime (t) for the normalized lifetime distribution P(t) is defined as 

POD 

(t) = dt P(t) t. (42) 







The random (exponential) lifetime distribution F = ke extremizes the information 
entropy 

/■oo 

S^ = - dt P(t) log[P(t)] (43) 
Jo 

where k = {t)~^ (see, for example, |124j ). One measure of the extent to which a 
calculated lifetime distribution P(t) characterized by mean lifetime (t) deviates from 
exponentiality is the entropy deficit |124] 

ASp = ^Random " ^(P) (44a) 

POD 

= 1 + log(t) + dt P(t) log[P(t)]. (44b) 
Jo 

(c) Flux through DS 



The flux through the DS is obtained by integrating the flux form da (cf. eq. (31 )) over 
the DS. As the flux form is simply the phase space volume element associated with 
the 'activated complex', the reactive flux (p can be computed by uniform (random) 
sampling of the DS phase space (see Appendix [B| . 
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The associated reactive volume 20 x s is the total phase space volume on the energy 
shell % = traced out by all the trajectories passing through the dividing surface 
(in both directions, hence the factor of 2). Except for a set of measure zero, each 
trajectory returns to the DS, and only contributes to the reactive volume up until the 
gap time. 

(d) Reactant density of states 

The classical density of states p{E) of the reactant region at energy H = is obtained 
by calculating a discretized approximation to p{E) as a function of energy for < 0, 
fitting p{E) to a polynomial function in and evaluating the fitted p{E) at -E = 0. 
The discretized approximation to p{E) is obtained by randomly sampling phase points 
inside a suitably chosen hypercube, and binning the energies for sampled points with 
E < 0. (See Appendix |B|) 

Although in all cases the phase space volume (integrated density of states) N{E) is 



found to be finite as — )■ 0, for the 3 DoF isokinetic Hamiltonian (eq. (20)) it is not 
clear from our numerical results whether or not the density of states actually diverges 
as — )■ 0, and we have not been able to evaluate this limit analytically. For the 



isokinetic thermostat with n = 4 DoF (eq. (23)), our numerical results suggest that 
p{E) is finite as E ^ 0. 

In Appendix |X| we show analytically for isotropic harmonic potentials that, while 
p{E) diverges logarithmically as E ^ for n = 2 DoF, for n > 4 DoF, with n even, 
p{E = 0) is finite. 

The value of the density of states p{E) at = is of some interest, as equality between 
the reactive density of states determined as the product of the mean gap time with the 
reactive flux and the total reactant density of states p{E = 0) is a necessary condition 
for ergodicity of the thermostat dynamics. 

(e) Thermostat dynamics 

It of course important to examine the effectiveness of the Hamiltonian isokinetic ther- 
mostats defined here as thermostats; that is, to assess how well time-averaged co- 
ordinate distributions evaluated over a single thermostat trajectory reproduce the 
Boltzmann distributions associated with the relevant temperature parameter f3. 
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We therefore pick an initial condition at the coordinate origin on the DS with random 
momenta and H = 0, and propagate it for a long time [t = 20000, many periods of the 
harmonic oscillator modes). Trajectories are integrated in Mathematica |125j using 
the function NDSolve with the SymplecticPartitionedRungeKutta method and fixed 
step size, At = 0.01. By binning coordinate values over such a long trajectory, we 
obtain a discretized probability distribution function that can be compared with the 
desired canonical (Boltzmann) coordinate distribution. 

Moments of powers x'^ and y'' obtained as time averages over the trajectory can also 
be compared with thermal averages at temperature T. 



B. Numerical results: 3 DoF 



1. Lifetime distributions 



Computed lifetime distributions for the 3 DoF thermostats H121, H321 and H521 are 
shown in Figure |2} It can be seen qualitatively from these plots of log[P(t)] versus t that, 
following initial transient decay, the lifetime distributions become more exponential as the 
temperature decreases (i.e., as /3 increases). The entropic measure of the deviation defined 



in eq. (44) was computed for (renormalized) lifetime decay curves with transients excluded 
(that is, we take t > (t)). The resulting values of AS are shown in Table |lj the AS values 
reflect the qualitative observation that the lifetime distributions become more exponential 
as (3 increases. 



2. Phase space volumes 

Numerical values for the flux, mean gap time and reactant phase space volumes are 
given in Table [ij The magnitude of the reactive flux (p decreases with temperature (as /3 
increases), as does the inverse gap time. Smaller inverse gap times therefore correlate with 
lifetime distributions that are more nearly exponential (cf. [SI])- 

Note that energy surface volumes for the 3 DoF systems are not shown in Table |Tj The 
reason for this omission is that, on the basis of our numerical computations (results not 
shown here), we cannot be sure whether or not p{E) diverges for 3 DoF systems as — )■ 
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Mean gap time 


Flux 


Reactive volume 


Energy surface volume 


AS 


H121 


16.57 


6.978 


231.28 




0.034 


H321 


48.88 


0.773 


75.61 




0.026 


H521 


130.41 


0.280 


72.99 




0.019 


J121 


12.69 


41.490 


1053.36 


1053.48 


0.037 


J321 


38.51 


1.536 


118.31 


118.66 


0.021 


J521 


101.60 


0.334 


67.87 


69.47 


0.010 



TABLE I: Computed mean gap times, fluxes, reactive phase space volumes, energy surface volumes 
and entropy deficits for lifetime distributions for 3 DoF and 4 DoF model Hamiltonians. Details 
of the computations are discussed in Appendix [Bj 

(cf. Appendix |A]). As we do not have an analytical proof that the density of states at E = 
is finite, we cannot rule out the possibility of a divergence. 

3. Thermostat coordinate distributions 

Coordinate distributions for the 3 DoF Hamiltonian isokinetic thermostat are shown in 
Figures |3} |4] and |5j The histogrammed distributions are obtained by time-averaging over 
a single trajectory, while the solid red curves are the associated canonical (Boltzmann) 
distributions. Moments {q'') for q = Xi are shown in Figure [gJ 

Qualitatively, it is apparent from these results that the effectiveness of the Hamiltonian 
isokinetic thermostat, as judged by the similarity of trajectory-based and Boltzmann co- 
ordinate distributions and moments, increases with temperature. Recall that the lifetime 
distributions become more nearly exponential as temperature decreases. 

Despite the possibility of an infinite energy shell volume, the Hamiltonian isokinetic 
thermostat for 3 DoF does in fact serve to thermalize the 2 uncoupled harmonic modes 
quite effectively. Note, however, that even for the long trajectories considered, numerical 
distributions obtained for the thermalizing (y) coordinate are not symmetric about y = 0. 
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C. Numerical results: 4 DoF 



1. Lifetime distributions 

Computed lifetime distributions for the 4 DoF thermostats J121, J321 and J521 are shown 
in Figure [7} Values of the lifetime distribution entropy deficit AS* are shown in Table [!} As 
for the 3 DoF systems, the 4 DoF lifetime distributions become more exponential as the 
temperature decreases. 



2. Phase space volumes 

Numerical values for the fiux, mean gap time and reactant phase space volumes for 4 
DoF systems are given in Table [1} Both the magnitude of the reactive fiux (p and the inverse 
gap time decrease with temperature (as (3 increases), while the lifetime decay curves become 
more nearly exponential as /3 increases. 

Comparison between reactive volumes and total energy surface volumes shows that both 
J121 and J321 systems satisfy (at least, within numerical error) the necessary condition for 
ergodicity, while the J521 system shows a minor deviation from equality. 



3. Thermostat coordinate distributions 



Distributions for cordinates Xi and y computed for the 4 DoF Hamiltonian isokinetic 
thermostat are shown in Figures |8] and [9| respectively. The histogrammed distributions 
are obtained by time-averaging over a single trajectory, while the solid red curves are the 
associated canonical (Boltzmann) distributions. Moments {q'^) for q = xi are shown in 



Figure 10 



As for 3 DoF, the effectiveness of the Hamiltonian isokinetic thermostat, as judged by 
the similarity of trajectory-based and Boltzmann coordinate distributions and moments for 
the harmonic oscillator coordinates, increases with temperature. Recall that the lifetime 
distributions become more nearly exponential as temperature decreases. 
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VII. SUMMARY AND CONCLUSIONS 



In this paper we have investigated the phase space structure and dynamics of a Hamilto- 
nian isokinetic thermostat. By design, ergodic thermostat trajectories at fixed (zero) energy 
generate a canonical distribution in configuration space [21 [Ml [731 El]- The physical poten- 
tials studied consist of a single bistable mode (the thermalizing degree of freedom [711 172] ) 
plus transverse harmonic modes. Although these modes are not coupled in the physical 
potential, the potential for the Hamiltonian thermostat is obtained by exponentiation of the 
physical potential, which introduces coupling between the modes. 

Interpreting the bistable mode as a reaction (isomerization) coordinate, we are able to 
establish connections with the theory of unimolecular reaction rates [IHl [HI I116[ I121j . in 
particular the formulation of isomerization rates in terms of gap times |26l [77] . In Thiele's 
general formulation ^77j, the gap time is the time taken for a reactive trajectory initiated on a 
dividing surface in phase space to return to the surface. (Such phase space dividing surfaces 
in multidimensional systems have been defined and computed using normal form theory [551 ^ 
[nS].) The distribution of gap times for a microcanonical ensemble initiated on the dividing 
surface is of great dynamical significance; an exponential distribution of the lifetimes for 
the reactive ensemble is usually taken to be an indicator of 'statistical' behavior. Moreover, 
comparison of the magnitude of the phase space volume swept out by reactive trajectories as 
they pass through the reactant region with the total phase space volume (classical density 
of states) for the reactant region provides a necessary condition for ergodic dynamics. If the 
total density of states is appreciably larger than the reactive volume, the system cannot be 
ergodic. 

We have computed gap times, associated lifetime distributions, mean gap times, reactive 
fiuxes, reactive volumes and total reactant phase space volumes for model systems with 
3 and 4 DoF. The symmetry of the model potentials studied means that in all cases the 
dividing surface is defined by a single condition on the thermalizing coordinate, y = 0. The 
thermostats were studied at three different temperatures. For 4 DoF, the necessary condition 
for ergodicity is approximately satisfied at all three temperatures. For both 3 and 4 DoF 
systems, nonexponential lifetime distributions are found at high temperatures (/3 = 1, where 
the potential barrier to isomerization is 1/2 in the same units), while at low temperatures 
(/3 = 5) the lifetime distribution is more nearly exponential. 
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We have quantified the degree of exponentiahty of the hfetime distribution by computing 
the information entropy deficit with respect to pure exponential decay. From the standpoint 
of unimolecular reaction rate theory, the decay becomes more "statistical" at lower T (smaller 
flux). This finding is in accord with the early observations of Deleon and Berne [5l] on 
isomerization dynamics in a model 2-mode system. 

We have examined the efficacy of the Hamiltonian isokinetic thermostat by computing 
coordinate distributions averaged over a single long trajectory initiated at random on the 
dividing surface. For the parameter values used here, coordinate distributions are more 
nearly canonical (Boltzmann-like) at lower temperatures. 

It remains for future work to establish more quantitative correlations between dynamical 
attributes of the Hamiltonian thermostat and thermostat effectiveness. 
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Appendix A: Phase space volume and classical density of states — analytical results 
for harmonic potentials 



An important question arising in numerical investigation of the dynamics and phase space 
structure of Hamiltonian isokinetic thermostats concerns the behavior of the density of states 
p{E) in the hmit E ^ 0. Specifically, it is important to know whether or not the density of 
states diverges in this limit. It can be difficult to establish convergence of p{E) to a finite 
value as -E — 7- on the basis of numerical calculations of the phase space volume N{E). 

In this Appendix, we evaluate analytically the phase space volume N{E) and the associ- 
ated density of states p{E) = dN{E)/dE for a Hamiltonian of the form 



n 



(Al) 



with 



Hi) 



Li=i 



(A2) 



and even dimensionality n. Physically, the potential (A2) corresponds to n uncoupled 



degenerate harmonic oscillators; we determine the phase space volume for the exponentiated 
harmonic potential appearing in the Hamiltonian ( |A1[ ). 

We start by evaluating the phase space volume N{E) enclosed by the energy shell Ti = 



E < 0. The potential (A2) is spherically symmetric, and so depends only on r, the radial 



coordinate in n-dimensional configuration space. Let r = r{E) be the value of r for which 
the exponentiated harmonic potential is equal to E, 



^/-2\og[-2E]. 



(A3) 



Note that, as — from below, r{E) — )■ oo. At fixed r, < r < r, the magnitude of the 
momentum Ivrl is 



1^1 = V2E + e-^V2 



V4^g{f2-r2)/2 _ I 



(A4a) 
(A4b) 



If V{x,n) is the volume of a ball of radius x in n dimensions. 



V{x, n) 



r[f + 1] 



(A5) 
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then the total phase space volume N{E) is given by the integral 



'•*'(^) dV(r n) 

N(E)= I dr V(|7fLn) V (A6a) 

dr 







In (A6b) 



r[| + i]2 " 

with 

rf{E) 

Xn = J dr nr''~\-'^'^l^ ^^-r-M^ _ ' _ ^^^^ 

For even dimension n, this integral can be evaluated explicitly. 
For n — 1 DoF, the integral X2 is 

X2 = 2-e-^"'/2(2 + ^2)_ (Ag) 

In terms of the energy we have 

N{E) = 27r2(l + E{2-2 log[-2E])), (A9) 

and the limiting value of the phase space volume is therefore finite 

\imN{E)\n=2^27r''. (AlO) 

Despite the fact that, when N{E) is plotted as a function of E, the derivative of N(E) 
appears to be finite at £■ = 0, the density of states p{E) is 

p{E)^^^^-4nHog[-2E] (All) 

which diverges as -E — )■ 0. This divergence is difficult to identify in a plot of N{E) versus 
E even when the form of N{E) is known, and is also difficult to see in plots of numerically 
determined p{E). 
For n — 4, 

X4 = e-*^' (r"^ + 6f^ - 16e^ + 14^ + 2. (A12) 

In terms of we have 

N{E) = (1 + 16E + 28^2 _ 24^2 iog[-2^] + SE^ log[-2^]2) , (A13) 
and the associated p{E) is 



p{E) = Stt^ (1 + 2£; - 2£;iog[-2£;] + £;iog[-2£;]^) 
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(A14) 



which is finite as — )■ 0. In fact, we have 



hm p(E) = 8n^ ~ 779.27. 



(A15) 



Figure 11 shows the ratio of the numerically determined p{E) for the 4 DoF harmonic 



potential to the theoretical expression (A14), over the energy range —0.1 < E < 0. The 



numerical p{E) is obtained by randomly sampling phase space points inside a suitably chosen 
hypercube, and binning points according to the value of the Hamiltonian. The resulting 
histogram provides a discretized approximation to p{E). 

The value of p{E = 0) for the 4D HO potential obtained using the procedure discussed 



in Appendix is 773.43, which is within 1% of the exact result (A15). 



To see the structure of these expressions in the general case for n = 2m, m integer, note 
that the integral 



2m 



f(E) 



dr 2mr2'^-^e-""^'/2 



g(f2-r2)/2 _ ^ 



(A16) 



is a sum of terms of the form (cj constant) 



f{E) 

J. = cje-^rn-j)fy2 I ^2m-lg-irV2^ J = 0, 1, . . . , m 



(A17) 

The only term that contributes to N{E) in the limit — )■ (f — )• oo) is that with j = m, 
so that 

^2m9^ pco _ 2m / 9 \ 

(A18) 



lim N(E) 

E-^O 



.0 , dr r^rn-l^-mr^/2 ^ 

V[m + 1Y Jq V[m + l\\m 



The limiting value of the phase space volume is therefore finite for all positive even integers 
n = 2m. 

To find the general form for the density of states, p{E), it is necessary to evaluate 
dX2m,/d-E. Differentiating the integral, noting that the integrand vanishes for r = f{E), 
and using 



dE 

df 



-f2/2 



we have 



dX, 



2m 



dE 



f{E) 



dr (2m) r 



2^2m-l -(m-l)fV2 



3(f2-r2)/2 



m— 1 



(A19) 



(A20) 



It is clear that p{E) will diverge as — > only for the case rri = 1; for all other positive 
integer m, the derivative dN{E)/dE converges to a finite value as -E — )■ 0. 

We have not been able to derive corresponding analytical expressions for p{E) for the 
isotropic harmonic case for odd dimensions. 
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Appendix B: Phase space sampling and computation of dynamical quantities for 
isokinetic thermostats 



In this Appendix we provide additional details of the computations discussed in Section 



The potentials for the isokinetic thermostat Hamiltonians discussed in Sec. |III| are ob- 
tained by exponentiating the physical potential, which has the form of a one dimensional 
double well plus uncoupled harmonic modes. For comparison with the analytical results for 
the density of states for an exponentiated isotropic harmonic potential given in Appendix 
[A] we have also investigated the corresponding Hamiltonians numerically for 2,3 and 4 DoF. 



As mentioned in Section |VI[ the symmetry of our Hamiltonians means that the dividing 
surface (DS) between 'reactant' and 'product' is simply defined to be the symmetry plane 
y = for the double well potential. 

The results reported in this paper were obtained using algorithms coded in Python |126] ; 
calculation efficiency was not a prime objective and could no doubt be improved significantly. 
Many of these calculations required large numbers of samples and were therefore subject to 
a trade-off between sample size and accuracy. The size of production runs was estimated on 
the basis of results obtained from shorter trial runs. 

At the value of the Hamiltonian "H = 0, the energy surface extends to infinity in coordinate 
(q) space; the energy of the saddle point is negative. As the exponentiated potential is 
bounded from below, momentum components are bounded but coordinate components are 
unbounded. For phase space points with larger coordinate component magnitudes, the 
relative probability of having energy E < decreases. A phase space sampling region 
("box") is therefore used with appropriate ranges of the momentum components and an 
upper limit on coordinate components chosen to be large enough that the (small) neglected 
volume associated with points having E < does not affect the accuracy of the result to 
within a chosen tolerance at the sampling density used. The largest values for coordinate 
component magnitudes were taken to be either 4.0 or 5.0 (in the units used here). 

Phase space points are randomly sampled within the box, and the required energy surface 
volume calculated by estimating the fraction of points which satisfy the appropriate energy 
condition. Increased accuracy requires both a large number of sample points and a large box; 
there is clearly a trade off between calculation accuracy, box size (hence neglected volume 
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outside the box), sample size and processing time. We have aimed for an accuracy of one 
percent while maintaining as far as possible a common approach across all the parameter 
sets studied. A count was taken of points in the periphery (outer 10% of coordinates) of the 
box to check on the validity of this process. 
Details of the computations now follow: 

(a) Mean gap time 

The mean gap time is calculated as an appropriately weighted average over the dividing 
surface. We sample the phase space manifold defined by constraints y = 0, Py = 0. 
Other coordinate and momentum components {x,Px) were randomly sampled within 
a chosen box. The value of the Hamiltonian was calculated at each point, and sample 
points with energies E > were rejected. Points with E < were used to seed a 
trajectory from the dividing surface. The value of the saddle point momentum py was 
calculated (with py > 0) so that the total energy H = 0, giving a trajectory initial 
condition on the dividing surface. Note that the momentum py is only calculated after 
a point on the DS has been accepted. The trajectory was then integrated using a 
4th-order Runge-Kutta algorithm until it recrossed the surface y = 0. The variation 
of the energy along the integrated trajectory was monitored to check stability of the 
integration procedure. The maximum cutoff time for the trajectory was either 2000 or 
5000 time units (for different parameter sets). The gap time for each trajectory was 
recorded and the average calculated over ~ 100,000 trajectories. 

(b) Flux 

The directional flux is the volume of the dividing surface y = at = (cf. eq. 



(31 )). It is calculated by uniformly sampling the restricted phase space region as above 
and again calculating the value of the Hamiltonian to decide whether each sampled 
point is within the dividing surface (that is, has energy < 0). The fraction of points 
accepted times the volume of the sampling region yields the volume of the DS for 

n = o. 

(c) Reactive volume 

The reactive volume is obtained by calculating (twice) the product of the mean gap 
time s and the one-way flux 0. Since only positive values of the saddle point momentum 
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are used, only trajectories on the positive side of the dividing surface are used to 
calculate the mean gap time. For symmetric systems such as those treated here, an 
additional factor of 2 is all that is required to calculate the total reactive volume on 
both sides of the dividing surface. 

(d) Energy surface volume and density of states 

We wish to calculate the density of states p{E) at energy E = 0, where p{E) is the 
derivative with respect to energy of N{E), the total phase space volume enclosed with 
the energy shell "H = 0. In order to do so we determine the volume of a set of energy 
shells of thickness AE. This data can then be used to obtain a polynomial fit to p{E) 
and hence its value at E = 0. 

A phase space sampling region (box) of full phase space dimensionality is used, so 
that we sample in two more dimensions than the computations described above. We 
sample phase space points z = (tt, g) G with Hlz) < 0, and consider the region of 
phase space occupied by points with energies —k x AE < E < —{k — 1) x AE, with 
AE typically 10~^ and k a positive integer. The volume of each of these energy shells 
is estimated by determining the number of points inside it as a fraction of those inside 
the full box; the corresponding density of states p{Ek) at E^ = —{k — ^) x AE can 
then be obtained by dividing the volume of the k-th shell by AE 

The density of states p{E) can then be calculated using either of two approaches. The 
total volume of a set of /cmax such shells can be calculated, and the resulting energy 
volume N{E) for — fcmax x '^E < E < fitted to, for example, a polynomial in E. The 
derivative of the fitting function can then be evaluated at AE = 0. An alternative 
procedure is to fit directly the estimated values of p{E) to a polynomial in E and use 
the resulting fit to obtain p{E) at E = 0. Both methods should in principle give the 
same result. 

In the calculations reported here, p{E) was calculated for a set of 2 x 10^ shells 
with AE = 10~^, and a fifth-order polynomial fit to p{E) calculated over the range 



—0.02 < E <0. Noise in the fitted data (see Figure 12) can be reduced by smoothing 
(concatenating several energy intervals AE), and we have verified that the values of 
p{E = 0) we obtain are robust with respect to such smoothing. 
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Figure captions 



FIG. 1: (a) Section (3:2 = 0) through the physical potential <&(xi,X2,y), eq. (19), wi = 1, uj2 = \/2, 



a = 2. (b) Section (x2 = 0) through the exponentiated potential, — ^ exp[— 2/3<I>], (3 = 1. 

FIG. 2: Lifetime distributions for the 3 DoF Hamiltonian isokinetic thermostat. The lifetime 
distribution is derived from the distribution of gap times obtained by initiating trajectories on the 
incoming DS and propagating them until they cross the outgoing DS. (a) H121, /3 = 1. (b) H321, 
/3 = 3. (c) H521, /? = 5. 

FIG. 3: Coordinate distributions for 3 DoF Hamiltonian isokinetic thermostat, obtained by av- 
eraging over a single trajectory. Numerical distributions for the xi coordinate (histograms) are 
compared with the Boltzmann distribution (solid line), (a) H121, j3 = 1. (b) H321, (3 = 3. (c) 
H521, (3 = 5. 

FIG. 4: Coordinate distributions for 3 DoF Hamiltonian isokinetic thermostat, obtained by av- 
eraging over a single trajectory. Numerical distributions for the X2 coordinate (histograms) are 
compared with the Boltzmann distribution (solid line), (a) H121, (3 = 1. (b) H321, (3 = 3. (c) 
H521, (3 = 5. 

FIG. 5: Coordinate distributions for 3 DoF Hamiltonian isokinetic thermostat, obtained by av- 
eraging over a single trajectory. Numerical distributions for the y coordinate (histograms) are 
compared with the Boltzmann distribution (solid line), (a) H121, (3 = 1. (b) H321, (3 = 3. (c) 
H521, (3 = 5. 

FIG. 6: Moments of the distribution of the coordinate xi obtained using the 3 DoF Hamiltonian 
isokinetic thermostat (squares) are compared with those for the Boltzmann distribution (circles). 
Odd moments for the Boltzmann distribution are identically zero, (a) H121, (3 = 1. (b) H321, 
(3 = 3. (c) H521, (3 = 5. 
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FIG. 7: Lifetime distributions for the 4 DoF Hamiltonian isokinetic thermostat. The Ufetime 
distribution is derived from the distribution of gap times obtained by initiating trajectories on the 
incoming DS and propagating them until they cross the outgoing DS. (a) J121, (3 = 1. (b) J321, 
(3 = 3. (c) J521, (3 = 5. 



FIG. 8: Coordinate distributions for 4 DoF Hamiltonian isokinetic thermostat, obtained by av- 
eraging over a single trajectory. Numerical distributions for the xi coordinate (histograms) are 
compared with the Boltzmann distribution (solid line), (a) J121, (3 = 1. (b) J321, (3 = 3. (c) J521, 
(3 = 5. 



FIG. 9: Coordinate distributions for 4 DoF Hamiltonian isokinetic thermostat, obtained by av- 
eraging over a single trajectory. Numerical distributions for the y coordinate (histograms) are 
compared with the Boltzmann distribution (solid line), (a) J121, (3 = 1. (b) J321, (3 = 3. (c) J521, 
(3 = 5. 



FIG. 10: Moments of the distribution of the coordinate xi obtained using the 4 DoF Hamiltonian 
isokinetic thermostat (squares) are compared with those for the Boltzmann distribution (circles). 
Odd moments for the Boltzmann distribution are identically zero, (a) J121, (3 = 1. (b) J321, 
/3 = 3. (c) J521, (3 = 5. 



FIG. 11: Ratio of the numerically determined p{E) for the 4 DoF harmonic potential to the 



theoretical expression (A14), over the energy range —0.1 < E < 0. The phase space volume N{E) 
was determined by random sampling of a phase space hypercube using Np = 5 x 10^ points. 



FIG. 12: Number of points per energy bin versus energy E (width AE = 10^^, —0.02 < E < 0) 
for the J321 4 DoF Hamiltonian. The red curve is a fit to the numerical data using a 5-th order 
polynomial. The fit to the data yields a value p{E = 0) = 118.66. 
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